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We study the dynamics of a system made up of particles of two different species undergoing 
irreversible quadratic autocatalytic reactions: A -\- B 2 A. We especially focus on the reaction 
velocity and on the average time at which the system achieves its inert state. By means of both 
analytical and numerical methods, we are also able to highlight the role of topology in the temporal 
evolution of the system. 



I. INTRODUCTION 

The interest in systems undergoing reaction-diffusion processes is experiencing a rapid growth, due to their intrinsic 
relevance in an extraordinary broad range of fields [1]. 

In particular, a great deal of experimental and theoretical work has been devoted to the study of reaction-diffusion 
processes embedded in restricted geometries. This expression refers to two, possibly concurrent, situations: i. low 
dimensionality and ii. small spatial extent. 

In the first case, the spectral dimension d characterizing the diffusive behaviour of the reactants on the substrate 
is low (1 < d < 2), and the substrate underlying the diffusion-reaction lacks spatial homogeneity. This situation is 
able to model media whose properties are not translationally invariant and where the reactants perform a "compact 
exploration" [2] . These kinds of structures can lead to a chemical behaviour significantly different from those occurring 
on substrates displaying a homogeneous spatial arrangement. Indeed, while in high dimensions a mean- field approach 
(based on classical rate equations) provides a good description, in low dimension local fluctuations are responsible for 
significant deviation from mean- field predictions [3]. 

There also exists a variety of experimental situations in which reaction-diffusion processes occur on spatial scales 
too small to allow an infinite volume treatment: in this case finite-size corrections to the asymptotic (infinite- volume) 
behaviour become predominant. 

Here, differently from previous works, we explicitly examine finite size systems, i.e. no thermodynamic limit is taken 
[4-6]. All the quantities we calculate are hence finite, and we seek their dependence on the finite parameters of the 
system (volume of the substrate and concentration of the reactants). In particular, we study the dynamics of a system 
made up of two species particles undergoing irreversible quadratic autocatalytic reactions A-\- B ^ 2 A. All particles 
move randomly and react upon encounter with probability 1, i.e. the reaction is strictly local and deterministic. 
Notice that, allowing all the particles to diffuse makes the problem under study a genuine multiparticle-diffusion 
problem. The latter is generally quite difficult to manage due to the fact that the effects of each single particle do not 
combine linearly, even in the non-interacting case. For this reason the analytic treatment often relies on simplifying 
assumptions which, nevertheless, preserve the main generic features of the problem. In the past, autocatalytic 
reactions have been extensively analyzed on Euclidean structures [7], within a continuous picture attained by the 
Fisher equation [8,9] which describes the system in terms of front propagation. Evidently, this picture is not suitable 
for low-density systems, where front propagation cannot be defined. In order to describe also the high-dilution regime, 
here a different approach is introduced which, as we will see, works as well for inhomogeneous structures. This way, 
we are also able to highlight the role of topology in the temporal evolution of the system. 

In the following, we shall examine the concentration Pa(^) of A particles present in the system at time t and its 
fluctuations; from Pa(^) it is then possible to derive an estimate for the reaction velocity. Furthermore, we consider 
the average time r (also called "Final Time") at which the system achieves its inert state, i.e. A^^ = A". As we will 
show, r depends on the number of particles A" and on the volume V of the underlying structure. More precisely, for 
small concentrations of the reactants, we find, both numerically and analytically, that the r factorizes into two terms 
depending on A" and V, respectively. 

One of the most interesting applications of the Final Time is analytic [10,11]: as we show, r sensitively depends 
on the initial amount of reactant A' and, on low dimensional substrates (d < 2\ by reducing the dimension the 
sensitivity can be further improved. 
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II. THE MODEL 

We consider a system made up of TV particles of two different chemical species A and diffusing and reacting 
on a discrete substrate with no excluded volume effects. At time A^a(^) and A^b(^) represent the number of A 
and B particles, respectively, with N = Na + Nb- Being V the substrate volume, we define pA{t) = NA{t)/V and 
PB{t) = NB{t)/V as the concentrations of the two species at time t. 

Different species particles residing at time step t, on the same node or on nearest-neighbour nodes react according 
to the mechanism A ^ B ^ 2A with reaction probability set equal to one. Notice that the previous scheme also 
includes possible additional products (other than 2 A) made up of some inert species of no consequences to the overall 
kinetics. The initial condition at time t = is Na{^) = 1 (the Source), Nb{0) = N — 1, with all particles distributed 
randomly throughout the substrate. As a consequence of the chemical reaction defined above, NA(t) is a monotonic 
function of t and, due to the finiteness of the system, it finally reaches value N; at that stage the system is chemically 
inert. The average time at which NA{t) = N is called "Final Time" and denoted by r. 

The Final Time r is of great experimental importance since it represents the average time when the system is inert 
and therefore it provides an estimate of the time when reaction-induced effects (such as side-reactions or photoemission) 
vanish [12]. In this perspective, deviations from the theoretical prediction of r are, as well, noteworthy: they could 
reveal the existence of competitive reactions or explain how the process is affected by external radiation. 

Finally, notice that the autocatalytic reaction can also be used as a model for spreading phenomena: A{B) particles 
may stand for (irreversibly) sick (healthy) or informed (unaware) agents, respectively. For these systems a knowledge 
of the infection rate or information diffusion is of great importance [4,5]. 

III. AVERAGE FINAL TIME 

As previously said, r generally depends on the total number of agents N and on the size of the lattice y, while 
its functional form is affected by the topology of the lattice itself. The analytical treatment is carried out in the two 
limit regimes of high and low density. 

A. High-density regime 

When p = N/V 1, the substrate topology does not qualitatively affects results. We can assume that the set 
of A particles covers a connected region of the substrate whose volume expands with a constant velocity (depending 
on the density p and dimension d). In this case (and exactly in the limit p oo) the process can be described as 
the deterministic propagation of a wave front decoupled from the random motion of the agents. If we suppose the 
Source to be at the center of the lattice at time t = 0, at each instant the wave front is the locus of points whose 
chemical distance from the center is 2t -h 1. The connected region spanned by the wave front is entirely occupied by 
A particles, while B particles fill the remaining of the lattice. In particular, for a (i- dimensional regular substrate, the 
region where A particles concentrate is a c?- dimensional polyhedron [4,5]. 

In general, for a finite system, the average Final Time is r = /max/2, where /max is the chemical distance of the 
most distant point on the lattice, starting from the Source. On Euclidean geometries this yields r = I//4 for d=l and 
T = L/2 ioT d > 2. On the other hand, on inhomogeneous structures, the dependence on L is not so simple, since it 
involves taking the average with respect to all possible starting points for the Source. 

B. Low-density regime 

In the case of low density (p <C 1) the time an A particle walks before meeting a B particle becomes very large, so 
that the process is diffusion-limited. We adopt a mean-field-like approximation by assuming that the time elapsing 
between a reaction and the successive one is long enough that the spatial distribution of reactants can be considered 
random. In other words, the particles between each event have the time to redistribute randomly on the lattice and 
we neglect correlations between their spatial positions. Another consequence of the low concentration of reactants, 
is that we can just focus on two-body interactions since the event of three or more particles interacting together 
is unlikely. Notice that the high-dilution assumption, by itself, generally does not allow to apply the classical rate 
equations: when diffusion is involved also the substrate topology has to be taken into account. For this reason, in the 
following we will treat high and low dimensional structures separately. 
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High- dimensional structures {d > 2) Let us consider a given configuration of the system where A^^ and Nb particles 
are present. The probabihty for a given B particle to encounter and react with any A particle is just the trapping 
probability Ptrap{PA,t) for a particle, out of TV^, in the presence of Na traps, both species diffusing. Under the 
assumptions specified above, for high- dimensional substrates [1]: 

Ptr.p{pA,t) = Xd PAe-^'P^\ (1) 

where A^^ is a constant depending on the given substrate. Form the previous equation we can calculate the average 
trapping time for a B particle as Ttrap(pA) ^ P~a • 

Let us now introduce an early-time {t < rtrap(pA)) approximation for the trapping probability: PtrapiPA^t) ~ pNa, 
where p ^ is the probability that, after each reaction, two given particles first encounter at a given time (in 
general, this probability depends not only on the volume of the underlying structure, but also on the history of the 
system). This simple form for Ptrap(pA, t) allows us to go on straightforwardly. In fact, the process can be meant as an 
absorbing Markov chain, with N states (labeled with the total number of A particles: 1,2,..., A/"), and one absorbing 
state {Na = N); the chain starts from state 1. The transition matrix P can be written: the transition probability 
from a state A: to a state m as a function of N and p is: 



k,m 



[1 -(!-#] [{i-p)"] (2) 



m 



k 



for any TV and p. From P we can take the submatrix Q, obtained subtracting the last row and column (those 
pertaining to the absorbing state), and compute the fundamental matrix F = (1 — Q)~^. Now, by expanding to first 
order in p^ a, direct calculation shows that F is an upper triangular matrix given by 

Fk,m = {^^^^^ l^'^ (3) 

y \j k < m. 

The mean time r required to reach the absorbing state N, starting from state 1 is given by the sum of the first row 
of F: 

r{N,V)='-y^ (4) 

^^'''•'^ p m{N -m) N^oo N ^ ' 

where 7 = 0.577... is the Euler-Mascheroni constant. The last result is in perfect agreement with numerical simulations 
and also emphasizes r factorization. 

Low- dimensional structures {d < 2). For low dimensional structures the dependence on TV found above is not 
correct. The reason is that a non-linear cooperative behaviour among particles emerges. 

Let us define {tn) as the average time elapsing between the (n — l)-th first encounter among different particles and 
the n-th one. This time just corresponds to the average time during which there are just NA(t) = n particles in the 
system. In our approximation {tn) is proportional to the trapping time rtrap{n/V) in the presence of n mobile traps 
diffusing throughout a volume V [6]. For compact exploration of the space {d < 2), Ttrap(^/^) ~ (V/n)^/^. This 
result was derived for infinite lattices, nonetheless, it provides a good approximation also for finite lattices, provided 
that the time to encounter is not too large. From rtrap we obtain {tn) as the average trapping time of the first out of 
N — n = Nb particles, that, for rare events, is just 

W^V^'^^. (5) 

with logarithmic corrections in the case d = 2. The time r can therefore be written as a sum over n(n= 1,2,..., AT— 1) 
of {tn)- Now, by adopting a continuous approximation, we obtain for d < 2 [6]: 



t{N, V) « V^/^ 



d 



(2 - d)N 



+ N-^'^ (logiV + if2/d) + 0{N-^) 



(6) 



where Hm is the harmonic number. In particular, the leading-order contribution for a one-dimensional system {d = 
d=l)\s 

r{N,V)^^. (7) 
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For a two-dimensional lattice {d = d = 2) 



r{N,V) ^ VlogV 



log TV + 7 
TV ' 



(8) 



Notice that the factorization in Eq.(|6]) is consistent with Eq.(|4j): in both cases, the factor containing the dependence 
on V represents the average time for two particles to meet. 
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FIG. 1: Scaling of r with the linear size of the system for a one-dimensional chain (blue circles), a Sierpinski gasket (black 
triangles), a T- fractal (red squares), and a three-dimensional cubic lattice (green diamonds) on a double- logarithmic scale. The 
number of reactants is fixed at A/" = 1024 for all systems. The spectral dimension for the Sierpinski and the T-graph is d ^1.365 
and d ~ 1.226, respectively. Dotted lines highlight the low-concentration regime (L ^ 1), corresponding to a power law for all 
systems. For the one-dimensional chain, the linear high-concentration regime is also pointed up. 



As can be evinced from Fig. 2, for small densities all the data collapse; moreover, in that region, the fit coefficients 
introduced are in good agreement with theoretical predictions. 




FIG. 2: Rescaled Final time versus number of particles N for the Sierpinski gasket (left) and the two-dimensional lattice (right). 
Different symbols and colours distinguish different sizes, as explained by the legend. The line provides the best fit in very good 
agreement with Eqs. (|6| and (|8]), apart from sub-leading corrections in the (marginal) case d — 2. 



For low densities, the standard deviation (Jr displays a dependence on N and L analogous to r; for high densities, 
cFr becomes vanishingly small, in fact the process becomes deterministic. 

As anticipated in Section 1, experimental measures of r are useful in monitoring trace reactants [6]. In the high- 
dilution regime, our results show that r = /j(TV)^jC^) and therefore, once the substrate size is fixed, the initial 
amount of reactant can be expressed as A/" = f^^(j / g^iy)). 
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A proper estimate of the sensitivity of this method is provided by the derivative dN/dr: the smaher the derivative 
and the larger the sensitivity. As can be evinced from Fig. 3, which displays numerical results for N and dN/dr ^ the 
smaller the concentration and the better the sensitivity of this technique. This makes such technique very suitable for 
the determination of ultratrace amounts of reactants, which is of great experimental importance [13]. Interestingly, 
dN/dr also depends on the substrate topology: when d < 2 and at fixed V, the sensitivity can be further improved 
by lowering the substrates dimension. Conversely, when d > 2^ dN/dr ceases to depend on d. 




FIG. 3: Log-log scale plot of the reactant amount N (left panel) and its derivative dN/dr (right panel) vs Final Time r. As 
shown in the legend, different substrate topologies (with approximately the same volume) are compared. Lines are guide to 
the eyes. 



IV. TEMPORAL EVOLUTION 

In this section we deal with quantities depending explicitly on time t. First of all, we consider the concentration 
pA of A particles present at time t. Due to the irreversibility of the reaction taken into account, pA is a monotonic 
increasing function; more precisely it is described by a sigmoidal law, typical of autocatalytic phenomena [7]. 

As shown in Fig. 4 the curves NA{t) grow faster, and saturate earlier, with increasing d {N and V being fixed). This 
is consistent with the meaning of the spectral dimension d: it describes the long-range connectivity structure of the 
substrate and the long-time diffusive behaviour of a random walker on the substrate. More precisely, for (i < 2, the 
number of different sites visited by each walker grows faster as d increases, and analogously the number of meetings 
between walkers. 

For d > 2 (e.g., d = 3 in the figure), A^a(^) is independent of d and is fitted by a pure sigmoidal function. Also 
notice that deviations between curves relevant to different topologies are especially important at early-times, while at 
long times they all agree with the pure sigmoidal curve. This result is consistent with the existence of two temporal 
regimes concerning diffusion on low-dimensional structures [1]. As a result, the topology of the underlying structure 
is important only at early times, while, at long times, the system evolves as expected for high-dimensional structures. 

Within the analytic framework developed in the last section, it is possible to derive some insights into the temporal 
behaviour displayed by NA{t). Being t{n) the average time at which the number of A particles reaches value n, 
recalling Eq. ([5|) we can write 

From which A^a(^) = ^~^(A^a), whose numerical solution provides an S-shaped curve consistent with data obtained 
from simulations. 

As for transient lattices, the easy form obtained for Ptrap(PA, t) ^ pNa and the assumption of a uniform distribution 
for agents positions, allow to write a Master equation for the number of A particles in the system: 

NA{t + 1) = NA{t) + {N- NAimi - (1 - (9) 

To first order in p: NA{t + 1) = Cp{NA{t)), being Cp a logistic-like map, with a repelling fixed point in (/'(O) = 
1 + Np)^ and an attracting fixed point in N {f'{N) = 1 — Np). Since Np ~ p <C 1, the increment of NA{t) at each 
time step is very small (of order p), and we can take the evolution to be continuous. Thus we obtain 

PA{t) = e^P\e^P' + A/- - 1)-^ (10) 

which is in good agreement with numerical results (Fig. 4). 



6 



1 

0.9 
0.8 

PJ 
0.6 
I 0.5 
0.4 

0.2 
0.1 



♦ d=3,V=40^ 




. Sigmoidal Fit 


\ 1 i 1 


* d=2,V=256^ 




♦ d^=1.365,V=(3^V3)/2 




♦ d^=1.226, V=3^°+1 


/ / 
















f / : 








/ 1 * 






1 


** 1 * 


^ J 













10^ 



10^ 

t 



10° 



10" 



FIG. 4: Normalized number of A particles NA{t)/N vs time t for a system made up oi N — 128 particles embedded on different 
structures, as explained in the legend. The best fit for the cubic lattice a pure sigmoidal function (see Eq. (5)), shown by the 
green line. The latter also provides the best fit for the long time behaviour of NA{t)/N on low dimensional substrates. 
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FIG. 5: Reaction velocity Fluctuations a a and concentration pA versus time for a system of = 128 particles diffusing on 
a Sierpinski gasket; three different generations (depicted in different colours) are shown. Notice tv < ta- 



From pA{t) one can derive the rate of reaction v{t) = dtNA{t) which represents the reaction velocity. As you can 
see from Fig. 5, in agreement with the theoretical predictions, v{t) is an asymmetrical curve exhibiting a maximum 
at a time denoted as ty^ obviously corresponding to a flex in NA{t). Interestingly, ty scales with the volume of the 
structure according to ty ~ y^/d ^j^i^]^ ig ^he same dependence shown by r. Moreover, at ty the population of the 
two species are about the same {NA{ty) = Nsity) = N/2). 

Hence, the efficiency of the autocatalytic reaction is not constant in time but, provided the number N of particles 
is conserved, it is maximum when the number of B particles is about N/2. From Eq. we can derive a similar result 
for the variance aA{t) of the number of A particles present on the substrate. Interestingly, fluctuations aA{t) peak at 
a time t^ which, again, depends on the system size with the same law as r; notice that t^ > ty. 
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V. CONCLUSION 

We introduced an analytic approach to deal with autocatalytic diffusion-reaction processes, also able to take into 
account the role played by particles discreteness and substrate topology. Within such framework, we derived in the 
low-density regime, for both fractal and Euclidean substrates, the exact dependence on system parameters displayed 
by the average Final Time, also highlighting how topology affects it. In particular, the case d = 2 is marginal. Exact 
results are also found for Euclidean lattices in the limit of high density. 

Theoretical results concerning the average Final Time find important applications in analytical fields, where mea- 
sures of r are exploited for detecting trace reactants. Our results suggest that the sensitivity of such technique is 
affected not only by the reactant concentration, but also by the topology of the structure underlying diffusion. 
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